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ABSTRACT 

The damping on the fluctuation spectrum and the presence of thermal velocities as 
properties of warm dark matter particles like sterile neutrinos imprint a distinct sig¬ 
nature found from the structure formation mechanisms to the internal structures of 
halos. Using warm dark matter simulations we explore these effects on the structure 
formation for different particle energies and we find that the formation of structure 
is more complex than originally assumed, a combination of top-down collapse and 
hierarchical (bottom-up) clustering on multiple scales. The degree on which one sce¬ 
nario is more prominent with respect to the other depends globally on the energy of 
the particle and locally on the morphology and architecture of the analyzed region. 
The presence of shells and caustics in warm dark matter halos is another important 
effect seen in simulations. Furthermore we discuss the impact of thermal velocities on 
the structure formation from theoretical considerations as well as from the analysis 
of the simulations. We re-examine the assumptions considered when estimating the 
velocity dispersion for warm dark matter particles that have been adopted in pre¬ 
vious works for more than a decade and we give an independent estimation for the 
velocities. We identify some inconsistencies in previous published results. The relation 
between the warm dark matter particle mass and its corresponding velocity dispersion 
is strongly model dependent, hence the constraints on particle mass from simulation 
results are weak. Finally, we review the technical difficulties that arise in warm dark 
matter simulations along with possible improvements of the methods. 
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1 INTRODUCTION 

Independent studies and observations of both small and 
large scale structure are presently challenging the otherwise 
widely embraced CDM model. The so-called missing satel¬ 
lites problem (e.g. Klypin et al. 1999 Moore et al. 1999), 


where observations of galaxies do not map the abundance 
of substructures that are produced in CDM cosmologies is 
a serious drawback of the model. Furthermore, at smaller 
scales, the density profiles of galaxies show large cores (e.g. 


de Blok et al.||2001| |Salucci et al.||20T2| |Kuzio de Naray"fc] 

Kaufmann|2011| ) that have not been reproduced by the sim¬ 

ulations. The failure to replicate in CDM simulations pure 
bulgeless galaxies which are observed in an important frac¬ 
tion (Kormendy et al. 20l0) adds to the problem. 

While several recipes have been proposed in the attempt 
of ameliorating these issues (e.g. Navarro et al. 1996; Mar- 
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tizzi et ak1|2013| |Mashchenko et aI7||2QQ8[ |Pontzen V Gov-] 


ernato||2012| ), most of them introducing baryonic physics 

processes, current studies conclude that even including re¬ 
peated baryonic outflows, large cored galaxies are not found 
in the simulations ( Marinacci et al.]|2Q14 ), although this is 
still highly debated in the literature. 

The aforementioned situations, where the CDM model 
proves deficient in explaining the observations, are demand¬ 
ing further investigation. The WDM models, with sterile 
neutrinos leading as most probable particle candidates have 
been well studied and discussed in the literature in the past 
thirty years with an increased interest in the last few years 
(see the highlights of Daniel Chalonge workshops and Col¬ 
loquiums 2011-2013 for latest developments in the WDM 


field ( |de Vega &; Sanchez|2 01l[ de Vega, Falvella & Sanchez 


2012 )). 


It has been shown recently (Destri, de Vega V Sanchez 


2012) that modeling the quantum pressure of fermionic par¬ 
ticles ( Weinberg||1962 Muccione L Pfenniger||200'6 ) on the 
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other hand, one can reproduce the expected cores in dwarf 
galaxies, known to be dark matter dominated. 

Moreover, the recent detection (Bul bul et al.|2014 Bo¬ 
yarsky et al.| [2014) of a 3.55 keV unidentified emission line 
both on the data from XMM spectrum of galaxy clusters 
and Chandra can be a hint of sterile neutrino decay. 

Since particles in warm dark matter models have differ¬ 
ent intrinsic properties from the cold dark matter particle 
candidates, the effect of these particle on structure forma¬ 
tion and evolution is expected to be qualitatively different 
on both large and small scales. 

Notwithstanding the difficulties in modeling properly 


ity is strongly dependent on the specific model adopted. The 


widely used formula for generating velocities (Bode et al. 


the neutrino particle, several attempts (e.g. Colombi, Do- 


delson &; Widrow|1996||Bode et al.|2001[|Macci6 et al.|2012| 


Ka mada et al.||2C)T3 Schneider et al.||2013 ) have been con¬ 

ducted with a successful outcome in solving some of the 
cases where CDM comes to an impasse. While the methods 
of CDM simulations have been accurately improved over the 
last decade, the WDM simulations encounter a number of 
difficulties in accurately describing the effect of such parti¬ 
cles on both large and small scales. In addition to the resolu¬ 
tion limitations that are met in the CDM case as well, WDM 
particles like neutrinos, for example, have a phase space den¬ 
sity limit, a Fermi-Dirac distribution and a thermal velocity 
dispersion. Referring merely to sterile neutrinos, these par¬ 
ticles decouple whilst relativistic. 

The effects of the initial velocities of the warm dark 
matter particle are expected to manifest themselves on small 
scale structure formation. The free streaming exponentially 
dampens the power spectrum of density fluctuations such 
that very few structures are formed below the damping scale. 
Conservation of the fine grained phase space density is ex¬ 
pected to set a maximum density that cannot be exceeded 
during the formation of structures with collisionless parti¬ 
cles. For a fermionic WDM particle, we can crudely define 
the coarse-grained phase-space density Q = p/a 3 , where p is 
the density and a is the velocity dispersion. This definition 
is only a good approximation for locally isotropic velocities 
where the density of particles is not strongly varying. 

Different numerical approaches have been used to ad¬ 
dress the impact of warm dark matter particles thermal ve¬ 
locities. Since the numerical resolution is strongly limited 
with respect to the physics, one knows that the phase space 
distribution sampling is anyway poor in space as well as in 
velocity space. The best compromise is to imprint the phys¬ 
ical particle velocity to the simulation particle, as common 
practice in galactic dynamics. The particle limited sampling 
is not a sufficient reason to entirely drop the velocity sam¬ 
pling by neglecting the thermal velocities as done in some 
previous works, while keeping nonetheless the power spec¬ 
trum cutoff implied by a non-zero thermal velocity (e.g. 
Schneide r et al.|2013 Governato et al.|2015 ). Nor is the fact 
that for some dark matter particles the thermal velocities 
are comparable or smaller than the bulk Zel’dovich veloci¬ 
ties. Even though it has been considered difficult to prescribe 
accurately initial thermal velocities in dark matter simula¬ 
tions, the importance of using them has been emphasized in 


previous studies like Colombi, Dodelson & Widrow (1996), 


Bode et al. (2001) and Melott (2007). 


2001) which sets such a relation is based on an assumption 
that overestimates the number of species at decoupling and 
in effect underestimates the value of thermal velocities. We 
will take the opportunity to discuss these assumptions and 
we will also provide a method for estimating thermal veloc¬ 
ities for fermionic, maxwellian and bosonic particles in both 
relativistic and non-relativistic regimes based on a different 
set of premises that takes quantum physics into considera¬ 
tion. 

Several analyses of warm dark matter simulations in 
the keV range conclude that the formation of structure is 
hierarchical, like in cold dark matter simulations. Tradition¬ 
ally, top-down structure formation means forming chrono¬ 
logically the biggest structures first and the smaller ones 
later, while bottom-up or hierarchical structure formation 
means the reverse scale order, making it difficult to describe 
a scenario in which both coexist. In fact it is well known 
since, e.g., Lin, Mestel, & Shu (1965) and Zel’dovich (1970) 


that typical structure formation proceeds in time first along 
pancakes, then filaments and then halos, mixing the large 
and small spatial scales at all times with different propor¬ 
tions. If we use this terminology in a broader sense, top- 
down describes the dominant long range effects on structure 
formation: sheets collapsing into filaments, collapsing into 
halos. Bottom-up hierarchical structure formation, on the 
other hand, describes dominant short scale effects, mergers 
of both early formed and later formed halos. We will examine 
how both of these mechanisms of structure formation show 
up in the warm dark matter simulations presented here. 

Additional constraints coming from peculiar features 
may be considered. In cold dark matter models, during the 
hierarchical evolution caustics are being wrapped inside ear¬ 
lier generations of the merging history, making them invisi¬ 
ble in some cases even at high resolutions. However, |Cooper| 
et al. (2010) show using cold dark matter simulations that 


In the absence of a tested universal mechanism of pro¬ 
duction for the warm dark matter particle, the relation be¬ 
tween the particle mass and its corresponding thermal veloc¬ 


accretion mechanisms of stars and dark matter clumps and 
disruption of the latter can produce concentering shells that 
resemble those observed in NGC 7600. In the warm dark 
matter simulations, as we will see, the shells and caustics 
are more visible, especially at high redshift, where the top- 
down formation occurs. 

Constraints on the mass of a warm dark matter parti¬ 
cle from Lyman-a Forest, cosmic weak lensing, gamma-ray 
bursts, etc. (e.g. Markovic &; Viel|2014 de Souza et al.|20l3 ) 
give a lower limit in the few keV range. To study the effects 
of warm dark matter on structure formation, we have, how¬ 
ever, explored a larger mass interval, focusing on the region 
where these effects are more prominent while fairly balanced 
by the resolution. 

The paper is structured as follows: in Section 2 we ex¬ 
plain the theoretical reasons for using the thermal compo¬ 
nent of velocities in the simulations. Subsection 2.2.1 shows 
how the common used formula for generating velocities in 
warm dark matter simulations ( Bode et al.||2001 ) is conjec¬ 
tured from hypothetical assumptions. Section 2.2.2 presents 
a different approach in estimating thermal velocities from 
the particle mass. In Section 2.3 we discuss some of the in¬ 
consistencies found in previous studies. Section 3 describes 
the parameters used in our simulations while Section 4 shows 
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the results found from analyzing the simulations. At last, we 
present our conclusions in Section 5. 


2 INITIAL CONDITIONS OF 

COSMOLOGICAL SIMULATIONS 

2.1 Velocities in the initial conditions of 
cosmological simulations 

The initial conditions of most CDM and WDM cosmologi¬ 
cal simulations have often initial thermal velocity taken as 
strictly zero, with the argument that at the finite initial 
redshift the thermal velocity of CDM particles is small in 
regard of the bulk flow following from Zel’dovich’s prescrip¬ 
tion. We argue below that this practice is numerically incon¬ 
sistent with the actual problem of describing a collisionless 
fluid of finite phase space density /. For structure forma¬ 
tion, the distribution of the dark matter particles in velocity 
space is most important, as stressed by Colombi, Dodelson 
k Widrow (1996). 


Indeed, integrating Newton’s equation of motion for a 
set of particles in a force field g is equivalent to solving in 
a Lagrangian way with discrete mass particles the collision¬ 
less Boltzmann equation with the characteristics method. In 
conventional notations the Eulerian description of the phase 
space volume conservation reads, 


df df Of n 

— T v • — -1- sr • — _0 

dt <9x <9v 


(i) 


where g is the force field. The mass density p is the projec¬ 
tion of the phase space density on velocity space: 


t>,)= / 


P(M) = / d 3 u/(x, v,t). 


( 2 ) 


The mass density p generates the force field g by New¬ 
ton’s gravity. In a cosmological setup the mean density po 
is subtracted, 


g = 


GId 3 x' \p(x!,t) -po{t)\ ^_^ [3 
= G J A'(i 3 » [/(x',v,()-/ 0 (v,t)] 


( 3 ) 


So in this context using vanishing small thermal velocity 
already poses a consistency problem since / is of the form 
£(v — vo(x))p(x). This implies representing the system with 
a diverging / in a vanishing fraction of the phase space vol¬ 
ume, in other words, mass belongs only to an infinitesimally 
thin 3D sheet in 6D phase space. As / is conserved along 
a characteristics, it implies that this singularity must per¬ 
sist at all subsequent times. Methods to conserve arbitrar¬ 


ily high phase space density have been set up (Abel 2012 
Hahn k Angulo 2015), but this is not necessarily sufficient. 


Ideally a physically sound solution fo to this Boltzmann- 
Poisson system should not be numerically sensitive to the 
initial condition discretization. In practice it is known that 
the gravitational iV-body problem is exponentially sensitive 
to perturbations ( Miller| |l964), so the best that can be ex¬ 
pected in such simulations is that over an ensemble of sim¬ 
ulations with identical statistical initial conditions, results 
follow a reproducible statistical distribution. Detailed evolu¬ 
tion of particles is sensitive to perturbations, but the average 
evolution of an ensemble of particles is predictable. 


When /o is finite and differentiable everywhere, in other 
words when fo mathematically exists, the sound situation 
that should be used, a slight variation of /o, a fluctuation, 
will also follow the same set of equations, so, writing / = 
fo + fi , g = go + gi, where f± and gi are the differences 
between the reference fo and the perturbed solution /, and 
using the fact that fo is a solution of the system, we obtain 
the exact equations for the differences f ±, and gi: 

dt + 9x +So dv Sl dv Sl dv } 


gl = 


G f d 3 x'd 3 v fi(x, v, t) 


(x-x') 


( 5 ) 


So we see that f± follows the same left-hand side equa¬ 
tion than fo in the unperturbed field go, except that now 
the right-hand side contains a source term whose first term 
gi • ^ is the product of the force fluctuations gi times the 
gradient of the original fo in velocity space. Thus a vanish¬ 
ing zero thermal velocity for a set of particles supposed to 
represent a physical fo is not only suspicious since it corre¬ 
sponds to a delta function in velocity space, but also because 
the gradient is at least as singular as fo- In other words 
a zero initial thermal width is inconsistent with the initial 
assumptions, and susceptible to arbitrary strong amplifica¬ 
tion of perturbations, since the variational equations con¬ 
tain a diverging source term to first order when the initial 
thermal velocity is small. The only possibility to cancel this 
diverging term is either to have vanishing force fluctuations 
gi, which is exceptional when fi is non-zero, or that gi is 
orthogonal to which is also exceptional. The second or¬ 
der term gi • on the right-hand side can cancel the first 
term only when fi is proportional and opposite to fo , which 
is also exceptional. In summary dealing with diverging fo 
is inconsistent with the implicit assumption of regularity of 
the mathematical problem. 

It is instructive to compare how simulating collisionless 
fluids is differently approached in the fields of stellar and 
galactic dynamics. While in cosmology the physical colli¬ 
sionless fluids is assumed to consist of elementary particles, 
in galactic dynamics the fluid is composed of stars. In both 
cases the numerical simulation particles are order of magni¬ 
tudes more massive than the physical particles. In the CDM 
context a simulation particle velocity is seen as representing 
the bulk flow of a large ensemble of CDM particles, explain¬ 
ing why zero velocity dispersion has been assumed. In galac¬ 
tic dynamics in contrast it is well known that doing so would 
immediately cause huge gravitational instabilities, and that 
the correct way to perform collisionless galaxy simulations is 
to ascribe to the simulation particles the same velocity dis¬ 
persion as the stars. This is also required for respecting the 
virial theorem. While the use of velocities has been shown 
(Mel ott|2007 ) to be of importance for the CDM simulations, 
in WDM it becomes even more relevant (Colombi, Dodelson 
k Widrow||1996 ) since the particles do have intrinsic non¬ 
zero velocity dispersion. Collisionless fluids are particle mass 
agnostic, when the particle mass starts to be important is 
when the 2-body relaxation time is shorter than the system 
age. Recalling the Chandrasekhar 2-body relaxation time in 


an arbitrarily large homogeneous medium (Chandrasekhar 
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1942l|Henon|T973l l, 


Frel 


87rG 2 mp\nN ’ 


(6) 


where v is the velocity dispersion of N particles of mass 
m with average density p, it is obvious that this relaxation 
time is proportional to v 3 , 1/m, and 1/lniV, so is exactly 
zero when v = 0. Thus CDM simulations with initial zero 
v are technically collisional until numerical noise heats the 
particles to larger v. 

From a completely different side of physics, the maximal 
phase space density constraint set by Heisenberg’s inequal¬ 
ity, 


AxAp x ^ 


h 

2 


( 7 ) 


gives a minimum particle velocity dispersion, taking Ax — 
n -1 / 3 and m the particle mass, 


A ^™ 1/3 
-- 

2 m 


0.003 



1/3 

km/s . (8) 


Taking a velocity dispersion lower than this value violates 
Heisenberg inequality, while taking it only slightly larger 
means that the particles behaviour is governed by quantum 
physics, not classical mechanics as assumed in all cosmolog¬ 
ical iV-body simulations. 

In summary, adopting even a slight non-zero velocity 
dispersion in cosmological simulation is certainly a safer and 
more correct physical assumption than taking strictly cold 
initial conditions associated with a mathematically singular 
and inconsistent state leading to situations not under con¬ 
trol on the numerical viewpoint due to singular phase space 
density. 


2.2 Thermal velocities of warm dark matter 
particles 

Depending on the WDM particle physics and properties, 
different scenarios can be considered regarding the particle 
velocity dispersion as function of redshift. In many scenarios, 
particles are copiously created in ultra-relativistic conditions 
and very good thermal equilibrium. As the Universe expands 
they may subsequently decouple on the thermal point of 
view because their collisional relaxation rate becomes lower 
than the expansion rate, while still interacting with the rest 
of matter by gravitational coupling. If the particles do not 
decay their comoving number density is conserved, and if 
they follow the collisionless Boltzmann equation their phase 
space density is conserved too. But this later assumption is 
more fragile because some residual elastic collisional relax¬ 
ation processes can still decrease the effective phase space 
density by coarse graining. 


2.2.1 Inspection of a commonly used result 

A frequently cited derivation of WDM particle velocities as 
function of mass and redshift can be found in IBode et al.l 
(2001). In their Appendix A they recall that for a thermal 
relict particle X that decouples when relativistic, the abun¬ 
dance n x relative to photons is: 


nx = f43/4\ \ gx 

y Qdec ) y H / 2 


( 9 ) 


where gdec is the number of relativistic species present at de¬ 
coupling, and gx is the number of spin states of the particle. 
Connecting then the particle mass density px — mxnx with 
the cosmological parameters Qx = p x /Pc and h, where the 
critical universal density p c = 3H 2 /8 ttG , and the Hubble 
constant H = 100/ikms -1 kpc -1 , they derive, 


O h 2 ~ 115 gx vfi x 
X ~ g dec 1.5 key ' 


( 10 ) 


We confirm this equation when using n 7 = 413 cm -3 . Then 
the authors proceed to derive a velocity formula. Since the 
distribution function of fermions without chemical poten¬ 
tial p is proportional to (exp(e x /kTx) + 1) _1 , they point 
out that if the particles decouple from photons when still 
relativistic ex = ( Px c2 + m 2 c 4 ) 1 ^ 2 — me 2 can be replaced 
by pxc where px is the particle momentum. To keep phase 
space density constant clearly in the relativistic regime px 
must stay proportional to Tx- But obviously as the regime 
passes to non-relativistic this argument does not hold, the 
exact formula valid at all Tx is 


Px 


kTx 

c 


2 

+ 2 kTx - 


(ii) 


Therefore it seems incorrect to assume that px is propor¬ 
tional to Tx also at low Tx- The exact scaling from Eq. CD 
becomes p 2 x oc 2 kTx, or 0.5 mxv 2 x oc kTx, that is, the ki¬ 
netic energy ex, not the momentum, scales as temperature 
at all redshifts. 

Another problem is the derived constant for velocity. 
They assume that the distribution function scales as the 
non-thermal distribution (exp(i>/i>o) + l) -1 , and give with¬ 
out detail vq (in their Eq. (A3)), 


vo(z) 
1 + z 




/keYV 
\m x ) 


kms 


( 12 ) 

where z is the redshift. But eliminating Qxh 2 in this previ¬ 
ous equation using Eq. ( |10| ) (their Eq. (A2)), we obtain for a 
mx = 1 keV particle (rounding also to 2 significant digits), 


Vp(z) 
1 + Z 


0.12 



keV 

-kms 

m x 


-l 


(13) 


Thus we find gdec = 1000 (gx/1.5) 1 / 3 . This is too high a 
value for gdec to be endorsed, as mentioned by the authors, 
by large entropy producing processes. Since the value for 
gdec varies linearly with the mass of the particle in the given 
cosmological model (Eq. 10), it allows the elimination of gd ec 
from the final expression for velocity This high value used 
for gdec leads to a significant decrease in the particle veloci¬ 
ties, as shown in Table 1. 

In the minimal standard model the number of the full 
set of particles is ~ 107 while in the minimal supersymmet¬ 
ric standard model, the value is increased to ~ 229 (Pier- 


paoli et al.|1998 ). Previous studies like Colombi, Dodelson & 


1 The authors cite a value of 688 for the number of relativistic 

degrees of freedom at the time of decoupling of a 1 keV particle, 
while then using the value of 1000. In Viel et al. 2005 this latter 
value is used, although a rigorous calculation gives 903 as the 
exact value. 






























Widrow ( 1996| p use a value of ~ 100 for right-handed neu¬ 
trinos decoupling before the electroweak phase transition at 


very high temperatures, while Pierpaoli et ah (1998) assume 
a conservative reference value of 150 for both gravitino and 
a standard warm dark matter candidate like the massive 
neutrino. 

The lower value for the velocity adopted by jBode et al.| 


(2001) has been used in most WDM simulations thereafter. 


This value for #dec, however, is valid for a 1 keV particle only 
if we assume that dark matter is made entirely by these type 


of particles, as shown in Eq. (10). For the cases in which a 


certain warm dark matter particle represents only a fraction 
of the total dark matter content, this value is different and 


Eq. (12) needs to be scaled accordingly. This aspect has been 


overlooked in some simulation studies of mixed dark matter, 
providing misleading results as we will show in Section 2.3. 


The next line following Bode et al. (2001) Eq. (A3) 
states: “The rms velocity is 3.571uo”. Recalculating the rms 


— (f,v/vo 


+ 


velocity of the adopted distribution function / = (e 
l) -1 we find a slightly larger factor: 

/ 2 f^°4:7rv 4 fdv £( 5 ) 2 n2 . 

= ^ = ( 14 ) 

where f is Riemann’s function. The slight discrepancy (the 
9 digit) appears thus as a misprint. 

In Appendix A we find that the largest correction factor 
for the rms speed of a Fermi-Dirac distribution with respect 
to a Maxwellian distribution is 1.07, not ~ 3.6 as stated in 


Maccio et al. (2012). The difference comes entirely from the 


very non-thermal distribution. 


2.2.2 Another scenario for quantum semi-degenerate 
particles 


In the previous Bode et al. (2001) scenario, WDM parti¬ 


cles are treated as localized mass bullets following Boltz¬ 
mann’s equation. However, at creation time they are also 
assumed to be ultra-relativistic and in thermal equilibrium 
with radiation and the rest of matter, typically following a 
Fermi-Dirac distribution since the known stable particles are 
fermions. This entails that their quantum nature does play a 
role at birth, they are at least semi-degenerate. Phase space 
density is high enough for the distinction between classical 
and quantum particles to play a role. The non-local Pauli 
principle applies then, each particle “knows” about the state 
of each other. Now if phase space density is approximately 
conserved then particles remain semi-degenerate at all times, 
which is inconsistent with the usual assumptions applied at 
low redshifts that they behave as classical particles. 

The known neutrinos offer a good example that parti¬ 
cles are quantum objects instead of localized mass objects. 
Real neutrinos are in addition of being fermions also in a 
superposition of three mass states. Since mass states prop¬ 
agate at different velocities, with time relict neutrinos are 
actually in a superposition of entangled mass states increas¬ 
ingly spread apart. How gravitational interaction with mat¬ 
ter structures can destroy the coherence of these entangled 


2 Interestingly this is the paper cited by [Bode et al. (2001) as 
reference for production mechanisms of WDM and their relation 
to cosmology 
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states is a question that will need to be addressed in future 
works. 

Here we develop a procedure to calculate precisely 
the particle velocity valid in all relativistic regimes for 
fermions or bosons. The full distribution Fermi-Dirac or 
Bose-Einstein distribution reads (e.g., Padmanabhan 2002]) 

(2nh) 3 exp((e — pi)/kT) d= 1 ’ ^ ^ 

where p is the particle momentum, g the spin-degeneracy 
factor (of order 1 or 2), pi the chemical potential, e = 
Vp^ + m 2 c 4 — me 2 the particle energy, m the particle 
mass, and T the particle temperature. The comoving num¬ 
ber density is calibrated according to a neutrino-like sce¬ 
nario where the particles are once coupled to the photon 
background and in thermal contact, at a time where gravi¬ 
tational perturbations are still linear. 

First, the assumption that the chemical potential pi is 
constant and negligible is not necessarily valid for identi¬ 
cal fermions which are created in a half-degenerate state. 
The Pauli principle has for effect that identical fermions, 
even with negligible interaction (like the weak nuclear force 
for neutrinos), possess an effective exchange potential , also 
sometimes called exchange-correlation potential (e.g., Atkins 
& Friedman 2005). In quantum chemistry and Density Func¬ 
tional Theory (DFT) the exchange potential is well known to 
be essential in the Hamiltonian describing electrons around 
a nucleus, or electrons in materials, although the exact form 
in different contexts is sometimes not well known. In the cos¬ 
mological context the exchange potential changes the chem¬ 
ical potential as the spatial density of identical fermions 
changes. This effective repulsive interaction for fermions 
makes the collisionless assumption of free fermions much 


less obvious. In Pfenniger & Muccione (2006) the effective 


interaction of free fermions was illustrated by solving ex¬ 
actly the time-dependent Schrodinger equation for two free 
but identical fermions in 3D space starting as Gaussian wave 
packets. In the quantum regime (high phase space density) 
these wave packets effectively interact and are scattered due 
to the repulsive exchange potential. In the classical regime 
(low phase space density) the wave-packets follow, as ex¬ 
pected, a straight trajectory. 

In quantum statistical mechanics (H uang|198 7) the ex¬ 
change potential between two particles has a well known 
form dependent on temperature and distance r 


4 >( r ) = 


—kT log ( 1 =F exp ( —mkT— 


= -kT log ( 1 =F exp ( -27r-^ 


(16) 

(17) 


where the — sign applies for fermions and + for bosons, and 
A is de Broglie wavelength. For semi-degenerate particles A is 
of order of n -1//3 , so in a semi-classical description fermions 
“feel” a rapidly varying repulsive force from neighbouring 
particles, while bosons an attractive force. The reality of 
the exchange potential can be invoked to cast a doubt that 
the commonly assumed collisionless approximation for semi¬ 
degenerate particles is valid in the cosmological context. In¬ 
stead one should expect a local thermalization of identical 
particles on a short time-scale. 

To calculate the chemical potential evolution in the cos¬ 
mological context, one needs therefore an additional assump- 
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tion, besides number conservation. The particle momentum 
and kinetic energy can not be assumed conserved due to 
the global gravitational interaction. A reasonable assump¬ 
tion ( Trodden &; Caroll|2004 ) is that the expanding medium 
proceeds adiabatically, at least as long as gravitational clus¬ 
tering is linear. This means that the particle specific entropy 
can be taken as a conserved quantity. 

The entropy S expressed as a function of other thermo¬ 
dynamical variables reads (e.g., Padmanabhan|2002 Vol. I, 
Eq. 5.73), 


S=-(E + PV-[iN) , 


(18) 


where E is the total thermal energy, P the pressure, V the 
volume, [i the chemical potential, and N the number of par¬ 
ticles. The specific entropy s = S/N divided by Boltzmann’s 
constant k is a pure number 


s(T, n) 

k 


e + P 
= iTP l —:- m 


1 

kT 


(19) 


where e = E/V is the specific energy density and n — N/V 
is the number density. 

The thermodynamical quantities for fermions and 
bosons at all regimes can be calculated accurately by eval¬ 
uating numerically the relativistic Fermi-Dirac or Bose- 
Einstein integrals for particle density n, energy density e, 
and pressure P (e.g., Padmanabh an|2 002, Vol. II, p. 216) as 
functions of temperature T and chemical potential fi [[] : 




exp((e — fi)/kT) =L 1 


( 20 ) 


e(T,n) = 


4t rg f° 

~hF Jo 


P 2 e 


exp((e — ji)/kT) ± 1 


dp , 


( 21 ) 


P(T,n) 


4 ng r°° p 2 1 c 2 p 2 

h 3 J 0 exp((e — p)/kT) =b 1 3 e + me 2 ^ ’ 


( 22 ) 

where g is the number of distinct particle states, and e = 
y/p 2 c 2 + m 2 c 4 — me 2 the particle energy. In the integrands 
the + sign is for fermions, the — sign for bosons. The con¬ 
served particle density n(T, fi) is related to universal expan¬ 
sion by the scale factor a = 1/(1 + z), thus 


n(T(z),n(z)) = n 0 (l + zf , (23) 


while the constant particle entropy gives 

= 4.20183245 , (24) 

k k 

For Fermi-Dirac particles the solution of this system 
for no = 115cm -3 , m — 1 keV, g — 1 in the non-relativistic 
regime is : 

^ = -1.6202, ^ =5.6186 • 10 12 . (25) 

For a graphical illustration of these functions behaviour, see 


3 This part follows closely the calculations made infpfenniger &| 
[M uccione| ( |2006| >, but correct a mistake where the used entropy 
expression was only valid in the ultra-relativistic regime, or when 

[i = 0. 


Fig. |Bl| in Appendix B. For a couple of dex around this so¬ 
lution the scaling with n, g and m for T and v goes with 
good approximation as follow 2 : 


T = 2.0654 • 10“ 


= 0.2226 


1\ 


2/3 


/ me 2 \ 

115 cm -3 g ) \keV ) 

1/3 /_ 2 \ -1 


1 


115 cm -3 g 


f me 2 
\keV 


kms 1 . (26) 


When the regime becomes relativistic this approximation is 
no longer accurate. One can solve the pair of equations (23) 
and (24) in any situation. 

In comparison, for Bose-Einstein particles the solution 
for the same parameters is: 


^ = -1.2451, ^ = 8.1348 • 10 12 . (27) 

Around this solution the scaling with n, g and m for T and 
v goes approximately as: 


T = 1.4265-10“ 


= 0.1768 


115 cm -3 g 


n 1 \ 2 ^ 3 / me 2 \ 1 
115 cm -3 g) \keVy 

fe) kms ' 


K 


1/3 


(28) 


If Maxwell-Bolzmann particles are used in simulations one 
can also calculate the solution, replacing the ±1 in integrals 
by zero, and taking s(oo, 0)/k = 4. The velocity coefficient is 
found to be 0.20592 km s -1 , intermediate between the Fermi- 
Dirac and Bose-Einstein cases. The correction of quantum 
statistics with respect to a Maxwellian distribution remains 
thus small, as demonstrated in Appendix A. 


2.3 Power spectrum of the warm dark matter 
simulations 

Since collisionless physics does not depend on the particle 
mass, the power spectrum must directly depend only on 
the velocity distribution of the particles, which results from 
the particle production mechanism. Colombi, Dodels on V| 
Widrow (1996) and Bode et al. (2001) also emphasize this 


point. 

To compute the transfer function for WDM models 
the fitting formula suggested by Bode, Turok and Ostriker 
(2001) gives: 


T '10= pcDM =[1 + W“1 


2 i^i —10 /u 


(29) 


where a , the scale of the break, is a function of the WDM pa¬ 
rameters, which are function of the velocity, while the index 
v is fixed. People prefer however, to use the mass dependence 


instead of the velocity, using Eq. (13) as a conversion. 


4 An astute reader might notice that for classical massive neutri¬ 
nos (0.01 < mc 2 /eV < 2) the found temperature is much lower 
than the commonly quoted temperature of 1.9 K. Actually the 
1.9 K value is valid for massless neutrinos only. The difference 
comes from the misleading use of temperature as an equivalent 
concept for energy and vice versa, while the neutrino rest mass 
energy does not contribute to thermal energy. The proper mean¬ 
ing of temperature is the quantity that would be measured, in 
the case of real neutrinos, by a cosmic sized thermometer able to 
thermalize with the neutrino background. 
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Viel et al. (2005) (see also Hansen et al. (2002)), using 


a Boltzmann code simulation, found that v — 1.12 is the 
best fit for k < 5 h Mpc -1 , and they obtained the following 
expression for a: 


a = 0.049 


t m x \-iii 

AM 0 ' 11 ! 

/ 7 \ 1-22 

(A) 

v 1 keV / 1 

1 0.25 J 

\ 0.7 / 


/i -1 Mpc. 


(30) 

In the case of warm dark matter particles, the stream¬ 
ing velocity supresses the matter power spectrum P(k) and 
the formation of structure, on scales smaller than their free- 
streaming scale. A rough estimation of the free-streaming 


scale is given by Bond et al. (1980): 


7 2 tt 

&FS — T- 

Afs 


(i®) (I) • < 3I > 


Depending on the model for the properties of a certain par¬ 
ticle, there can be different expressions for the damping of 
the power spectrum (Abazajian & Koushiap pas|2006 


e-g.. 


but for the purpose of our present work and for easier com¬ 
parison with previous studies we use the expression given in 


Eq. (30) with the corresponding thermal velocities. 


This approach used for cutting the power spectrum is 
only valid however, for a scenario in which the whole dark 
matter content is made up by one specific dark matter par¬ 
ticle of a certain velocity. 


2.4 Caveat Emptor 


In this section we would like to summarize the findings of 
previous sections and discuss some of their implications. We 
want to stress that the assumed particle production model 
and physics strongly impact on the ascribed particle mass, 
while the initial velocity distribution and its corresponding 
power spectrum is the only really important initial parame¬ 
ter influencing the simulation results. As far as the physics 
behind the origins of the dark matter particles is concerned, 
the assumptions found in the literature can widely differ. 
Previously on section 2.2.1. we showed how the |Bode| 


et al. (2001) result for estimating velocities for neutrino like 


particles is based on arguments like entropy production and 
negligible chemical potential. The expression for velocities in 


Eq. (12) is based on a dependence of the number of species 


on the mass of the particle, such as to preserve the equiva¬ 


lence in Eq. (10). For the models with cold plus warm dark 
matter, or models with different particle mass, the value of 
the number of species should be adjusted accordingly. Many 
papers that study mixed particles simulations have omit¬ 


ted this readjustment for velocities (e.g. Anderhalden et al. 


2012). Eq. (12) has been reduced by the fraction with which 


a certain particle contributes to the total density, therefore 
leading to inconsistencies like having for a certain mass of 
a particle, different thermal velocities, depending solely on 
that fraction. Moreover, since the power spectrum cutoff de¬ 
pends on the velocity of the particle (not the mass), studies 
that use the cutoff for a velocity, but a different thermal ve¬ 
locity, given by a different model of particle production are 
not consistent. These results, although used for constrain¬ 
ing the mass of a particle in terms of detection experiments, 
should not be considered as accurate. 

As an alternative, we provide a different energy-thermal 
velocity correspondence based on number conservation and a 


Table 1 . Correspondence between particle mass m and rms ve¬ 
locity dispersion in literature for 0.2, 1 and 3.5 keV. The first 
column shows the value originally given in Bode et al. (2001), the 
second column shows the value obtained using g^ ec = 150 ( |Pier-| 
|paoli et al. 1998|) in Eq. |l2| , and the third one, the value given 
by our derivation. 


Mass 

Bode et al. 

Eq. (Il)x3.571 

Pierpaoli et al. 

This work 
Eq. (25) 

ke V / c 2 

km/s 

km/s 

km/s 

0.2 

0.366 

0.4032 

1.113 

1.0 

0.0429 

0.0225 

0.223 

3.5 

0.00806 

0.0230 

0.0636 


non-entropy production while taking into account the quan¬ 
tum pressure, but assuming a thermalization caused by the 
exchange potential. Entropy conservation by particles in the 
hot Big Bang is invoked by many authors, such as |Padman-| 
abhan (2002) or Weinberg (2008). From Eq. (26) which esti¬ 


mates the thermal speed of WDM particles, independent of 
the cosmological parameters, we have the following velocity 
dependence with the redshift: 


1 + z 


= 0.2226 


n l\ 1//3 / me 2 \ 
115cm- 3 g) (kev) mS 


(32) 

The difference between our and Bode et al. (2001) estima¬ 
tions is showed in Table. 0 for 0.2, 1 and 3.5 keV respec¬ 
tively, at redshift zero. The |Bode et al. (2001) speed for 
1 keV fermions out of equilibrium, used in most WDM sim¬ 
ulations, is 5 times lower than the value derived here. In 
general these differences cumulate their effect if simulations 
are started at higher redshifts, and are crucial not only for 
phase space density studies, but also for structure formation. 

Our finding affects the results and conclusions of previ¬ 


ous papers which were using Eq. (12) to constrain the mass 


of sterile neutrinos. This extends even to papers which did 
not include thermal velocities. The power spectrum studies 
based on the velocity of the particle are subject to the same 
difference in the velocity estimation (see Section 2.3). Also, 
when comparing the thermal velocities to the Zeldovich ve¬ 
locities, these factors weaken correspondingly the reason for 
ignoring the thermal velocities, against all the arguments 
presented in Section 2.1. 

Since our aim here is to describe typical qualitative 
effects on structure formation present in a broader range 
of energies, we will refer to the particles in terms of their 
velocity dispersion instead of their mass. Indeed the ther¬ 
mal velocity of a particle as its decoupling temperature at 
a certain redshift depends on the specific physics of particle 
production. That influences the ascribed mass of that par¬ 
ticle. More complex analysis of the decoupling theories for 
a certain particle may give a slightly different dependence 
between the thermal velocity at a certain redshift and the 
particle mass. 


3 SIMULATIONS SETUP 

We conducted several suites of iV-body simulations. All 
simulations have been performed once with pkdgrav-2, 
a treecode written by Joachim Stadel and Thomas Quinn 
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Figure 1. Illustrative density map of structure formation regions at redshift zero, from left to right, in CDM, WDM (0.3km/s) and 
HDM (2.3km/s). For a similar illustration of a full box see Maccio et al (2012). These simulations have not been used in this paper. 


Table 2. Details of the simulations 


Label velocities z\ cutoff box size N softening 
km/s eV Mpc/h pc 

CDM no 1 40 300 s 50 

WDM1 no 200 40 300 s 50 

WDM2 36.6 200 40 300 3 50 

WDM3 no 1000 40 300 3 50 

WDM4 4.6 1000 40 300 s 50 

WDM5 3643 200 30 256^ 100 


(Stadel 2001 Q a nd then using Volker Springel’s Gadget-2 


(Springel 2005 The initial conditio ns are generat ed with 


Ed Bertschinger’s GRAFIC2 package (Bertschinger 200 T 


Although some differences have been spotted between the 
two different codes, those differences are not qualitatively 
important where structure formation is concerned, there can 
be minimally spotted at very small scales. 

The simulations we have performed cover a range of ve¬ 
locities from 0.01 km/s to 10 km/s (3.5 keV to 15 eV) at red- 
shift zero. For illustration purposes, in Fig.^we show generic 
density maps of structure formation regions in CDM, WDM 
and HDM simulations. Particles that have ~ 0.1 km/s veloc¬ 
ity dispersion are in a transient regime from a predominant 
top-down structure formation scenario to a hierarchical one, 
showing both these trends. We have chosen one such simu¬ 
lation and compared it to a simulation of a colder particle 
more favored by the observational constraints and with a 
cold dark matter simulation. For the warm dark matter the 
simulations have been prepared with initial power spectrum 
consistent with initial velocities, and, for comparison, the 
same initial power spectrum without initial velocities. 

The simulation parameters are summarised in Table [2] 
The starting redshift for the simulations is Zi — 100 in order 
to ensure a proper treatment of the non-linear growth of 
cosmic structures. 

The cosmological parameters used are given by the 


5 http://hpcforge.org/projects/pkdgrav2/ 

6 http://www.mpa-garching.mpg.de/gadget/ 

7 http://web.mit.edu/~edbert/ 


WMAP7 results: O\=0.721, Q m =0.279, fl 6 =0.0463, h = 0.7 
and ag = 0.821, 

We start with running large scale simulations in cosmo¬ 
logical box of 40 Mpc/h, engaging 300 3 dark matter particles 
and one 30 Mpc/h box with 256 3 particles. We then select 
a region where the top-down halo formation is predominant 
and re-simulate it with an eight times higher resolution. 


4 SIMULATIONS ANALYSIS AND RESULTS 

4.1 Structure formation in WDM simulations 

Free streaming causes a delay in the formation of structure 
in the warm dark matter simulations. This delay depends 
on the energy, hence the velocity of the particle. The higher 
the thermal velocity of the particle, the later the filaments 
will reach the collapse, making it impossible for structures 
to be formed by redshift zero in hot dark matter scenarios, 
as illustrated in Fig. [T] right panel. 

The fragmentation of the filaments is observed in all N- 
body simulations of warm dark matter where the collapse 
is stimulated by the noise in the particle distribution (|Bode 
et al. | [2QQT] |Gotz &, Sommer-Larsen|[2002[ |Wang &; White 
2007| ) . This is especially observable at the characteristic grid 
or glass spacing. Above the free streaming scale the mass 
function is flattened to a value that closely matches the lu¬ 
minosity function of galaxies (assuming mass traces light). 
The length and the lifetime of the filaments depend on the 
energy of the particle. The higher the velocity dispersion of 
the particle, the larger will be the filaments and the longer 
they will be preserved. These can reach 20 Mpc in a 40 Mpc 
box and survive until redshift zero in the case of velocities 
of few km/s and above. 

In Fig. [2] the difference between high density regions in 
a CDM simulation, versus WDM at redshift 2.3 is shown. 
The picture displays the 2D projection of the 3D density 
map of the full simulation box. One can see that due to the 
free streaming, particles are concentrated in large spatial 
structures, delimited from each other by large low density 
regions, or voids, as opposed to the crowded web present in 
the cold dark matter simulation. 

It is well known that in CDM models, smaller halos 
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Figure 2. Plot of high density regions at redshift 2.3 in a 40Mpc/h simulation box with CDM particles in the left panel and WDM 
particles (model WDM2) in the right one, showing major topological differences 


collapse first and merge hierarchically into larger systems, 


as it is obvious in all high resolution simulations (Diemand 
&; Moore||2010 2011 Stadel et al.||2009 ). Furthermore, one 


finds that less massive halos are more concentrated, perhaps 
reflecting the fact that the density of the universe is higher 
at earlier epochs, since the CDM particles have an infinite 
phase space density. 

On the other end of the velocities spectrum, for HDM 
models, the structure formation is essentially top-down up 
to redshift zero, as can be seen in the right panel of Fig. [l] 
with just large filaments collapsing into few halos. 

As stated in the introduction, in the case of warm dark 
matter we see from our simulations that the structure forma¬ 
tion is more complex, a hybrid mechanism where both long- 
range and short-range effects are present, from long distance 
to nearest neighbours, from top-down to bottom-up. 

The top-down trend predominant in the early epochs 
in warm dark matter simulations has been missed in pre¬ 
vious works since it is difficult to observe it while ana¬ 
lyzing the snapshots of the simulations. For particles with 
velocities smaller than a few km/s the top-down trend is 
hidden by the hierarchical growth that dominates at later 
times. We have been able to see this effect in our simu¬ 
lations while watching movies made with a sufficient large 
number of snapshots. We stress that only movies convey 
the complexity of these multiscale hierarchical processes. 
Several movies of WDM structure formation, filament col¬ 
lapse and halo formation from this study can be found on 
a youtube playlist (https://www.youtube.com/playlist? 
list=PLnGS4wkSt Jlaqi3M9hTDaUzuZ-vs-Qg6i J^J 

As the movies show, in WDM simulations, structures 
form in a qualitatively different way from CDM models. The 


8 All the HD movies are on youtube and can be watched in- 
divadually on this channel https://www.youtube.com/channel/ 
UCEmQi8hDNW2emqGn-urtvpg Remember to adjust your settings 
to HD quality. Links for direct download can be provided on de¬ 
mand. For a short description of the movies and movies snapshots, 
see Appendix C. 


hybrid structure formation is more complex than what the 
traditional top-down/bottom-up dichotomy can categorize, 
as discussed in the fntroduction. 


• During the early stages one sees the formation of well 
contoured filaments. How early is this stage depends on the 
particle velocity. In our WDM2 simulations this happens in 
the interval 13 > z > 8. 

• In the higher density regions, usually situated at the 
intersection of such filaments, the first halos are formed 
through gravitational collapse. These halos continue grow¬ 
ing into larger ones by accreting particles from the disrupted 
filaments (Fig. |Cl| ). 

• In medium density regions, halos show a hierarchical 
formation trend. Small halos collapse first and then merge 
into bigger halos (Fig. |C2| . 

• In less dense regions, the ones isolated by voids and 
which have a very slow evolution, we have observed filaments 
that collapse very late. The top down formed halo survives 
without any mergers until redshift zero (Fig. |C3| . 

• Finally there is the more complex scenario in which 
we observe large halos formed earlier which merge together 
forming a large cluster (Fig. C4). 

• The filamentary-like structure is preserved until redshift 
zero, with new filaments forming in the low density regions 
as late as redshift z ~ 4. 


Looking closer, we have analyzed four different regions 
in our simulations and displayed them in four different 
movies. The characteristics of these regions are summarized 
in Table OH 

Our conclusion from analyzing these simulations is that 
there is only one correlation, between the time of the first 
collapse and the density reached in a certain region, and 
that depends only on the network morphology and archi¬ 
tecture of that region. The first halos collapse in the region 
where the density becomes ~ 2 x 10 3 times larger that the 
average density and almost ~ 3 x 10 3 times larger than the 
lowest densities present in that region at that epoch. In the 
simulations with particle velocity of 0.36 km/s (that mimic 
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Figure 3. Halo formation at the intersection of filaments. A zoom in projection shows that shells and caustics are visible in the not yet 
virialized WDM halo. 



Figure 4. A thin slice of the WDM halo formed Top-Down on the left and of a CDM one on the right. Very different internal structure, 
with shells in caustics in the WDM halo being more apparent 


Table 3. The properties of four different regions of simulation 
WDM2 displayed in the movies 


Label 

size 

first collapse 

average density 

highest density 

- 

box 

z 

critical 

critical 

lu.avi 

1/4 

10.13 

0.264 

477 

ld.avi 

1/4 

9.45 

0.258 

481 

ru.avi 

1/4 

10.77 

0.268 

480 

rd.avi 

1/4 

9.78 

0.258 

474 


200 eV), the first collapse appears just after redshift 10 with 
the first halos forming until redshift 4, while in the simula¬ 
tions with 0.04 km/s (that mimic 1 keV), the first structures 
would have been already formed by redshift 10. The first ha¬ 
los form at the high density region at the intersection of the 


filaments and then continue accreting matter. If in a certain 
region there are many filaments collapsing, then the halos 
will merge into bigger ones. 


Due to the free streaming velocity of the particles, the 
network configuration and architecture of a certain region is 
rapidly changing. When the density becomes higher in more 
isolated regions, the collapse occurs even later, after redshift 
4 and some of those halos do not suffer mergers (Fig. C3), 
so there are halos at redshift zero that have formed via a 
top-down scenario and did not grow through hierarchical 
mergers . This is an interesting result, since the observations 
show that a large fraction of halos in the universe have not 
suffered any mergers until redshift zero. 


Why a certain region has more a top-down or bottom-up 
formation history depends only on the spatial distribution 
of the filaments in a certain simulation. 
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Figure 5. A zoom in simulations with the same cutoff in the power spectrum having corresponding thermal velocities - WDM2(left) 
and no thermal velocities WDM 1 (right) at z=0 


A single halo simulated with different velocities can be 
seen in the movie halo.avj^] The 7 x 10 12 M© halo forms top- 
down at the intersection of the filaments and has 18 million 
particles in its 7*200 radius. 

These high resolution runs are 8 3 times more resolved 
in mass than the initial ones: the dark matter particle mass 
is rrtd — 2.72 x 1O 5 M0, where each dark matter particle has 
a spline gravitation softening of 355 pc. 

Although the WDM halos on galactic scales contain few 
bound substructures, one can see shells and caustics inside 
the virialized region which arise from the coherent infall of 
material along filaments and from the smooth surrounding 
regions. As the resolution increases, the presence of shells 
and caustics becomes more apparent. The early top-down 
formation of a halo at the intersection filaments is shown in 
Fig.[3]along with a zoom in its central region. One can clearly 
see the shells and caustics wrapped inside the 18 million 
particle halo. A thin slice projection of the warm dark matter 
halo and a cold dark matter one, clearly illustrates in Fig. [4] 
how strikingly different is their inner structure. 

4.2 Impact of thermal velocities on structure 
formation 

As stressed in Section 2, the use of thermal velocities in 
warm dark matter simulations is crucial, even if their value 

9 The movie can be found on youtube https://www.youtube. 
com/watch?v=s_H4dS0P27I A zoom in the halo can be watched 
at https://www.youtube.com/watch?v=zqVi9SSWmXM 


is comparable with the Zeldovich velocities at a certain red- 
shift. 

In Fig. [5] we show the differences at redshift zero be¬ 
tween the structures emerging in a region in two similar sim¬ 
ulations, WDM1 (without thermal velocities), and WDM2 
(with thermal velocities). Both simulations have the same 
size, the same power spectrum cutoff and the same initial 
redshift. The structure formation and evolution in these two 
simulations is shown side by side in movie cosmoboxall.avf^) 

We can see that although the position of the big struc¬ 
tures is not affected, below Mpc scales there is a memory of 
the grid in the simulation without velocities that is smoothed 
out when adding thermal velocities, as expected. Some of the 
very small halos formed in the simulation without velocities 
cannot be found in the simulation where thermal velocities 
are included. The lack of small halos in WDM simulations 
with velocities is of course a crucial feature hinting to resolve 
the discrepancy between the CDM simulations predicting 
too many subhalos in galaxy-sized halos in comparison with 
the observed number of dwarf galaxies around large galax¬ 
ies. Indeed WDM simulations without velocities still suffer 
from the infinite phase space density problem. 

For comparison, we have performed a suit of simula¬ 
tions that start with a cold distribution of particles, no 
power spectrum cutoff, but have velocities corresponding 
to lkm/s, lOkm/s, 200km/s and 700km/s. Even the early 
structure formation is qualitatively different from the warm 


10 The movie can be watched here: https://www.youtube.com/ 
wat ch?v=5txGwBRNC1U 
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dark matter simulations, confirming that the top-down col¬ 
lapse is induced by the damping of the power spectrum at 
small scales and not the thermal velocities. 


4.3 Technical aspects in simulating WDM 


The resolution limit poses even a more stringent problem 
in warm dark matter simulations then in cold dark matter 
ones. Indeed, in order to properly analyze a region of the 
simulation, multiple refinements of that region with higher 
resolution particles are used. This implies tracking the par¬ 
ticle backwards, from redshift zero to the initial conditions. 
Due to the large streaming velocities, particles that end up 
in a virialized halo at redshift zero come from a larger re¬ 
gion than in the CDM simulations, making it more difficult 
to reach high resolution simulations in WDM. 

The heavier the effective mass of our simulation parti¬ 
cles, the more prominent is the 2-body relaxation effect in 
small clumps (Eq. This problem is more stringent in the 
case of cold dark matter simulations, where an initial zero 
velocity is used. In the case of warm dark matter, this scales 
with the velocity of the particle, giving a smaller relaxation 
time for a smaller velocity. This is why for simulations in 
the keV range, where the streaming velocity is smaller, the 
top-down formation history has been barely observed. 


As recently shown by Gao, Theuns &; Springel (2015), 
methods like ’FoF’ used in analyzing cold dark matter simu¬ 
lations are proved to be insufficient in analyzing warm dark 
matter halos. We confirm this statement, finding that the ar¬ 
tificial fragmentation occurring along the filaments results in 
a high number of small halos with less than ten particles. 


5 CONCLUSIONS AND DISCUSSION 

We have performed several N-body warm dark matter sim¬ 
ulations within a large range of velocity dispersion, for the 
purpose of pointing out the effects on the formation of struc¬ 
ture. We have then focused on a regime where the resolution 
is better balanced by the velocity distribution. Some of our 
findings are summarized below. 

• In warm dark matter models, as our dark matter only 
simulations show, the structure formation follows a hybrid 
scenario in which both top-down and bottom-up scenarios 
have a saying. 

• The early structure formation in this warm dark matter 
models is essentially top-down, with large halos forming in 
the highest density regions, tracked at the intersection of fil¬ 
aments. The second level of top down formation of structure 
is occurring along single isolated filaments. 

• The biggest earlier formed halos will accrete matter 
from the filaments, while in small densities regions the merg¬ 
ers of smaller halos will result in a larger halo. 

• Later on, and depending on the morphology of the re¬ 
gion in which these halos formed, meaning mainly the den¬ 
sity and the layout of the filaments, they merge into bigger 
halos creating a hierarchical build-up. 

• The warm dark matter halos, especially the ones that 
did not suffer big mergers, show obvious shells and caustics. 

• The warmer the dark matter the more pronounced is 
the top-down effect and the more delayed is the collapse. 


• Albeit the numerical limitations we encounter as far 
as our warm dark matter simulations are concerned, we 
can conclude that an early top-down structure formation 
trend would be seen even in dark matter simulations with 
v < 0.05 km/s. For colder particles, this effect is hidden 
and wiped out by following abundant mergers resulting in a 
redshift zero distribution that seems in agreement with the 
hierarchical formation scenario. 

• The number of small satellites, as previously found, is 
visibly reduced in the WDM simulations with respect to the 
CDM ones. 

For a warm dark matter particle, as supported by the 
arguments adduced in Section 2, the thermal component of 
the velocity is important for different theoretical and prac¬ 
tical considerations. The strong dependence of the mass- 
velocity relation on the actual particle production model 
makes it difficult to constrain certain properties of the dark 
matter particle, including its mass. The impact that a cer¬ 
tain velocity dispersion is having on the structure formation 
and evolution on both small and large scales, as seen in sim¬ 
ulations cannot be used as a strong constraint on the mass 
in the absence of a universal model for particle production. 
Furthermore, we have shown that there have been some in¬ 
consistencies in previous studies with respect to the use of 
velocities in the simulations, that lead to ambiguous results. 


The baryonic physics may play an important role in 
the actual formation and evolution of halos, hence the ne¬ 
cessity of further exploring these effects in simulations. High 
redshift observations of halos could be used in comparison 
with complex baryonic warm dark matter simulations in 
constraining the mass of warm dark matter particles based 
on their formation and merger history. 

The baryonic processes that we have not included in 
the simulations must play an important role in the structure 
formation. Previously Gao & Theuns (2007) show a crucial 
difference in the collapse of a filament that contains both 
gas and dark matter in a 3keV simulation, with respect 
to the cold dark matter case. In the WDM case, the stars 
form inside the filament, before the halo forms. This trend 
where stars form in filaments continues for 1.5 keV particles 
up to redshift z ~ 2 resulting in stringy “chain” galaxies 


that remain to be confirmed by observations (Gao, Theuns 
fc Springel|2015 ). 


The smoother space distribution in the warm dark mat¬ 
ter scenario may allow baryons to condense coherently in 
a smooth potential halo, providing favorable conditions for 
forming disk-like galaxies. However, a much higher resolu¬ 
tion that the one available in present simulations is needed 
to explore this effect. 
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Figure Al. The ratios between the velocity dispersions 
ctfd/omb and cbe/cmb with respect to the temperature 


APPENDIX A: VELOCITY DISPERSION 
DEPENDENCE ON TEMPERATURE IN 
FERMI-DIRAC AND BOSE-EINSTEIN 
DISTRIBUTIONS 

As mentioned in Sect. |2.2TT) we check the correction to veloc¬ 
ity dispersion that should be applied to a Maxell-Boltzmann 
distribution when the physical system follows a given quan¬ 
tum statistics. 

The energy of a particle as a function of momentum p, 
valid in all regimes (relativistic and non-relativistic), is 


e (p) — y/p 2 c 2 + m 2 c 4 — me 2 , 

(Al) 

and the velocity 


pc 

v(p) — = . 

y/m 2 C 2 + p 2 

(A2) 

The corresponding ID spherical distributions for 
Dirac, Maxwell-Boltzmann and Bose-Einstein cases 

Fermi- 

are: 

_ 47rp 2 

FD exp (e(p) / kT) + 1 

(A3) 

j. Airp 2 

/MB " exp (e(p)/kT) 

(A4) 

„ 4:irp 2 

/BE ” exp ( e(p)/kT ) - 1 

(A5) 

Integrating over all p, we obtain the normalization constant 

POO 

S= fdp. 

Jo 

(A6) 

For each case respectively the second moment are 


2 1 r 2 < * 

^FD = ^ J P JFB dp , 

(A7) 

2 1 2 f A 

%B = ^ J P JMB dp , 

(AS) 

2 1 r 2 , , 

C^BE = ^ J P /BE dp . 

(A9) 


Computing these integrals by numerical quadrature, we find 


the ratios between the velocity dispersions ctfd/ctmb and 
ctbe/ctmb with respect to temperature. The result is plotted 
in Fig. |Al| In any situation the Fermi-Dirac velocity disper¬ 
sion is not significantly different from Maxwell-Boltzmann’s, 
differing by at most ~ 6.5%, while the Bose-Einstein velocity 
dispersion differs more, up to ~ 27%. The highest differences 
occur at low temperature, corresponding to low redshifts. 
This is not such a dramatic correction as the factor 3.571 
invoked in Maccio et al. (2012), but can still be significant 
for high precision cosmology works. 


APPENDIX B: DETAILED DERIVATION OF 
THE RESULTS IN SECTION 2.2.2 


Using these expressions inserted into Eq. (19), the specific 
particle entropy becomes 


s 

k 



3 / 2 V2^+y(5g+4y) 

Z ~ 1 exp(?/)=l=l 
y 1 / 2 ^2q+y(q+y) 
Z ~ 1 exp(^)=Ll 


HZ) , 


(Bl) 


where Z = exp (p/kT), q = mc 2 /kT and y = 
(i/p 2 c 2 + m 2 c 4 — mc 2 )/kT. Thus s is a function of the re¬ 
duced dimensionless variables Z and q only, and not of g, m 
and physical constants explicitly. 

In the ultra-relativistic regime when the energy of parti¬ 
cles is comparable or higher than the rest mass energy, parti¬ 
cles and their antiparticles can be created in equal number, 
so any chemical potential should cancel to a high degree. 
Then s/k at fi = 0 becomes a constant. The closed form 
expressions are, 


lim 

T—> oo 


s(T, 0) 
k 


7 7T 4 

135 C(3) 


4.20183245, 


(B2) 


for fermions, and 

lim 0), _ 4^ _ 3 00157071, (B3) 

T—>oo k 45 C(3) 


for bosons, where £ is Riemann’s function, and £(3) ~ 
1.20205690. 

The particle velocity at all regimes can be derived 
from the relativistic particle kinetic energy e(T, pt) = 
e(2»/n(2» = m 2 c 4 — me 2 and that relativistic 

momentum is related to velocity by u 2 /c 2 = l/(l+m 2 c 2 /p 2 ). 
Eliminating p yields, noting Y = e/mc 2 , 


v 2 (T,p) = 1 _ Y (2 + Y) 

c 2 (1 + Y) 2 (1 + Y) 2 


(B4) 


The second form is numerically more precise at low energy. 
The non-relativistic and ultra-relativistic expansions read, 
respectively, 


% PS 2Y -3Y 2 + 4Y 3 - ... (B5) 

c 2 
2 

% PS 1 - Y~ 2 + 2 Y~ 3 - 3 Y~ 4 + ... (B6) 

c z 

As stated in Section 2.2.2, the conserved particle density 
n(T, fi) is related to universal expansion by the scale factor 
a — 1/(1 + z) and therefore 


n{T)z),ii{z)) = n 0 (l + z) 3 , (B7) 
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Figure Bl. Density n and entropy s/k as functions of temper¬ 
ature T and chemical potential fi as given in Eq. ( |B7| > and ( |B8| >. 
The light-grey surface is the entropy and the dark-grey ones are 
the density at two redshifts z = 10 9 on the left (relativistic), and 
z = 0 on the right (non-relativistic), while no = 115 cm -3 , g = 1, 
and m = IkeV/c 2 . The intersection of the level curves yields the 
solution of Eq. ( |B7| > and ( |B8| ). 


while the constant particle entropy gives 


s(T(z),n(z)) _ s(oo, 0) 


k 


k 


= 4.20183245 


(B8) 


For a given particle density no, redshift z and particle mass 


m the non-linear Eq. (B7) and (B8) can be solved with a 
non-linear equation solver for T and g. The functions are 
univalued and level curves of n and s intersect once, so any 


combination of T and g gives a single solution (see Fig. Bl). 
Actually the constant level curves n and s expressed with 
the variables logg and Z intersect almost at right angle: 
n(logg, Z) depends most rapidly on log q, and s(logg, Z) 
depends most rapidly on Z, so finding a solution for logg 
for n at constant Z and then a solution for Z at constant q 
for s, and repeating until satisfaction could be a method to 
find a solution. Since the thermodynamic functions involve 
integrals, a fast numerical integrator is handy, since sev¬ 
eral indefinite integrals must be evaluated at each iteration. 
To perform this we used Maple 18 which includes a non¬ 
linear multidimensional function root solver, and evaluate 
quickly numerical integrals with the NAG library algorithm 
DO1AMC0 When T and g are found for a given particle 
mass and degeneracy factor g , all the other quantities like 
v 2 can be derived by plugging these values in the functions, 
which may require again few numerical integral evaluations. 
The results are presented in Section 2.2.2. 


li 


The Maple script is available on request. 
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Table Cl. Description of the movies accompanying the paper 


Label 

cosmoboxvel.avi 

cosmoboxall.avi 

lu.avi ld.avi ru.avi rd.avi 

halo.avi 

halozoom.avi 


Description 

Movie of full-box WDM2 simulation 

WDM1 and WDM2 full-box simulations side-by-side showing the effect of thermal velocities 
A zoom in the 1/4 of the WDM2 simulation 
A 7 X 10 12 Mq 18 X 10 6 particles high-resolution refined halo from WDM5 
A zoom in the refined halo focused on the central region where the shells and caustics can be observed 



Figure Cl. A zoom in a region from the WDM2 simulation, showing the evolution of a halo which forms top-down at the intersection 
of the filaments and then starts accreting matter 



Figure C2. A zoom in a region from the WDM2 simulation showing how small halos formed later that merge hierarchically in a larger 
halo 



Figure C3. An early formed halo which doesn’t suffer mergers 



Figure C4. A large high density region with many filaments where the halos formed early on via top-down collapse are merging in a 
violent manner creating a larger cluster 
























